*******************************************************************************
*
* Discrete duration model + selection models
*
*******************************************************************************

eststo clear

*Beardsley's original robustness check for selection effects using a discrete duration model
use book-dyadyeardata, clear
eststo, title(Original): biprobit (newcrisis=manip2 manip2_t nomanip2 nomanip2_t prevcris2 viol2 crisdur2 jointdem victory2 contig2 prevcris2_t viol2_t crisdur2_t jointdem_t victory2_t contig2_t) (med2=viol2 ethniccomp2 crisdur2  lncapratio geostr2 medprev2) if _t<3650, cluster(dyadno) 

* --> Rho is estimated to be very small and insignificant = no sign for selection effect


*repeat the same selection model estimation with restricted model (Model 3 in Table 1)
* + allow for changing baseline hazard probability
eststo, title(Restricted + time polynomials): biprobit (newcrisis=manip2 manip2_t nomanip2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 prevcris2_t c._t c._t#c._t c._t#c._t#c._t) (med2=viol2 ethniccomp2 crisdur2  lncapratio geostr2 medprev2) if _t<3650, cluster(dyadno) 

* --> Rho is estimated to be very small and insignificant = no sign for selection effect

*estimate the model without selection equation
eststo, title(W/o selection equation): probit newcrisis manip2 manip2_t nomanip2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 prevcris2_t c._t c._t#c._t c._t#c._t#c._t if _t<3650, cluster(dyadno) 

* --> as expected, given insignificant and small rho: substantively the same coefficient estimates

esttab, unstack mtitle drop(/:)
esttab, star(* 0.1 ** 0.05 *** 0.01) b(3) se(3)  mtitle nodepvar  /*
*/ 		order(manip2 manip2_t nomanip2 nomanip2_t) noobs /*
*/		scalar("rho" "N" "ll") unstack  drop(/:) /*
*/ 		varlabels(manip2 "Manipulation"/*
*/ 		manip2_t "Manipulation * time"/*
*/		nomanip2 "Other mediation "/*
*/ 		nomanip2_t "Other mediation * time"/*
*/ 		prevcris2 "Previous crises"/*
*/ 		viol2 "Violence level"/*
*/ 		crisdur2 "Crisis duration"/*
*/ 		jointdem "Democratic dyad"/*
*/ 		victory2 "Victory"/*
*/ 		contig2 "Contiguity"/*
*/ 		prevcris2_t "Previous crises * time"/*
*/ 		viol2_t "Violence level * time"/*
*/ 		crisdur2_t "Crisis duration * time"/*
*/ 		jointdem_t "Democratic dyad * time"/*
*/ 		victory2_t "Victory * time"/*
*/ 		contig2_t "Contiguity * time"/*
*/ 		lncapratio "ln(Capability ratio)"/*
*/ 		ethniccomp2 "Ethnic component"/*
*/ 		geostr2 "Geostrategic salience"/*
*/ 		medprev2 "Previous mediation"/*
*/ 		_t "Analysis time"/*
*/ 		c._t#c._t "Analysis time²"/*
*/ 		c._t#c._t#c._t "Analysis time³"/*
*/		_cons Constant)
esttab using _table_A1.rtf, star(* 0.1 ** 0.05 *** 0.01)  b(3) se(3)  mtitle nodepvar  /*
*/ 		order(manip2 manip2_t nomanip2 nomanip2_t) noobs replace /*
*/		scalar("rho" "N" "ll") unstack  drop(/:) /*
*/ 		varlabels(manip2 "Manipulation"/*
*/ 		manip2_t "Manipulation * time"/*
*/		nomanip2 "Other mediation "/*
*/ 		nomanip2_t "Other mediation * time"/*
*/ 		prevcris2 "Previous crises"/*
*/ 		viol2 "Violence level"/*
*/ 		crisdur2 "Crisis duration"/*
*/ 		jointdem "Democratic dyad"/*
*/ 		victory2 "Victory"/*
*/ 		contig2 "Contiguity"/*
*/ 		prevcris2_t "Previous crises * time"/*
*/ 		viol2_t "Violence level * time"/*
*/ 		crisdur2_t "Crisis duration * time"/*
*/ 		jointdem_t "Democratic dyad * time"/*
*/ 		victory2_t "Victory * time"/*
*/ 		contig2_t "Contiguity * time"/*
*/ 		lncapratio "ln(Capability ratio)"/*
*/ 		ethniccomp2 "Ethnic component"/*
*/ 		geostr2 "Geostrategic salience"/*
*/ 		medprev2 "Previous mediation"/*
*/ 		_t "Analysis time"/*
*/ 		c._t#c._t "Analysis time²"/*
*/ 		c._t#c._t#c._t "Analysis time³"/*
*/		_cons Constant) nonotes  /*
*/ 		addnotes("cluster robust s.e., * p < 0.1, ** p < 0.05, *** p < 0.01") 



*******************************************************************************
* Graphical interpretation


*estimate the model with stata interaction to allow calculation with margins
* + allow for changing baseline hazard probability
probit newcrisis manip2 c.manip2#c._t nomanip2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 c.prevcris2#c._t c._t c._t#c._t c._t#c._t#c._t if _t<3650, cluster(dyadno) 

*we can calculate the predicted hazard probability of a new crisis
margins, at((mean) _all nomanip2=0 manip2=(0 1) _t=(0(100)3650))
marginsplot, xdim(_t) xtitle(Days since last crisis) ytitle(Hazard probability)  recast(line) recastci(rarea) ci2opts(fc(gs10%60) lw(none)) ci1opts(fc(black%60) lw(none)) plot1opts(lc(black) lp(dash)) plot2opts(lp(solid) lc(black)) legend(order(3 "No mediation" 4 "Manipulation")) title("")

* --> after ten years, the predicted probability of a new crisis is not even significantly higher
* --> no indication of a mediation dilemma

*calculate survival function for non-mediated cases
margins, at((mean) _all nomanip2=0 manip2=(0) _t=(0(365)3650))
matrix b0=r(b)		//save predictions as vector
matrix b0t=b0'		//transpose vector 
svmat b0t			//store vector as new variable
rename b0t1 hazardprobability_0		//rename
gen cumh0=sum(hazardprobability_0)	//calculate cumulative hazard probability
gen s0=exp(-cumh0)					//calculate survival probability

*calculate survival function for cases with manipulative mediation
margins, at((mean) _all nomanip2=0 manip2=(1) _t=(0(365)3650))
matrix b1=r(b)
matrix b1t=b1'
svmat b1t
rename b1t1 hazardprobability_1
gen cumh1=sum(hazardprobability_1)
gen s1=exp(-cumh1)

*generate time variable for plot
gen time_in_years=_n-1 in 1/11
gen time_in_days=time_in_years*365

*graph survival functions
twoway line s0 s1 time_in_days, ytitle(Survival probability) xtitle(Days since last crisis) legend(label(1 "No mediation") label(2 "Manipulation")) lp(dash solid) lc(black black)

* --> Substantively the same results as with Cox model
* ---> no indication of a mediation dilemma

